home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_400 / 440_01 / !best002.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-04-15  |  18.2 KB  |  732 lines
  1. /*==========================================================================
  2.  *
  3.  *  !BEST002.C                                       Tuesday, April 12, 1994
  4.  *
  5.  *  source file 2D and 3D matrix/vector library
  6.  *  supplementary source file 2 for The BESTLibrary
  7.  *
  8.  *  Authored by George Vanous with reference to the 2D and 3D vector library
  9.  *  by Andrew Glassner from "Graphics Gems", Academic Press, 1990
  10.  *
  11.  *==========================================================================*/
  12.  
  13.  
  14. /* ------------------------------------------------------------------------ */
  15. /* ----------------------------  INCLUDE FILES  --------------------------- */
  16.  
  17. #include <math.h>
  18. #include "!best002.h"                  /* include !BEST002.H in compilation */
  19.  
  20. /*==========================================================================
  21.  *                               3X3 MATRICES
  22.  *--------------------------------------------------------------------------*/
  23.  
  24. /*----------------------------------------------------------------------------
  25.  * C = A + B
  26.  * RETURNS: C
  27.  */
  28. matrix3 *m3_add(matrix3 *A, matrix3 *B, matrix3 *C)
  29. {
  30.   word i;
  31.   register word j;
  32.  
  33.   for (i = 0; i < 3; i++) {
  34.     for (j = 0; j < 3; j++)
  35.       C->element[i][j] = A->element[i][j] + B->element[i][j];
  36.   }
  37.   return( C );
  38. }
  39.  
  40. /*----------------------------------------------------------------------------
  41.  * RETURNS: copy of A
  42.  */
  43. matrix3 *m3_copy(matrix3 *A)
  44. {
  45.   matrix3 *B = NEWTYPE(matrix3);
  46.   register word i, j;
  47.  
  48.   for (i = 0; i < 3; i++) {
  49.     for (j = 0; j < 3; j++)
  50.       B->element[i][j] = A->element[i][j];
  51.   }
  52.   return( B );
  53. }
  54.  
  55. /*----------------------------------------------------------------------------
  56.  * C = A * B
  57.  * RETURNS: C
  58.  */
  59. matrix3 *m3_mul(matrix3 *A, matrix3 *B, matrix3 *C)
  60. {
  61.   word i;
  62.   register word j, k;
  63.  
  64.   for (i = 0; i < 3; i++) {
  65.     for (j = 0; j < 3; j++) {
  66.       C->element[i][j] = 0;
  67.       for (k = 0; k < 3; k++)
  68.         C->element[i][j] += A->element[i][k] * B->element[k][j];
  69.     }
  70.   }
  71.   return( C );
  72. }
  73.  
  74. /*----------------------------------------------------------------------------
  75.  * RETURNS: new 3x3 matrix initialized to a[3][3]
  76.  */
  77. matrix3 *m3_new(double *a[3][3])
  78. {
  79.   matrix3 *A = NEWTYPE(matrix3);
  80.   register word i, j;
  81.  
  82.   for (i = 0; i < 3; i++) {
  83.     for (j = 0; j < 3; j++)
  84.       A->element[i][j] = *a[i][j];
  85.   }
  86.   return( A );
  87. }
  88. /*----------------------------------------------------------------------------
  89.  * C = A - B
  90.  * RETURNS: C
  91.  */
  92. matrix3 *m3_sub(matrix3 *A, matrix3 *B, matrix3 *C)
  93. {
  94.   word i;
  95.   register word j;
  96.  
  97.   for (i = 0; i < 3; i++) {
  98.     for (j = 0; j < 3; j++)
  99.       C->element[i][j] = A->element[i][j] - B->element[i][j];
  100.   }
  101.   return( C );
  102. }
  103.  
  104. /*----------------------------------------------------------------------------
  105.  * B = transpose(A)
  106.  * RETURNS: B
  107.  */
  108. matrix3 *m3_transpose(matrix3 *A, matrix3 *B)
  109. {
  110.   register word i, j;
  111.  
  112.   for (i = 0; i < 3; i++) {
  113.     for (j = 0; j < 3; j++)
  114.       B->element[i][j] = A->element[j][i];
  115.   }
  116.   return( B );
  117. }
  118.  
  119. /*==========================================================================
  120.  *                               4x4 MATRICES
  121.  *--------------------------------------------------------------------------*/
  122.  
  123. /*----------------------------------------------------------------------------
  124.  * C = A + B
  125.  * RETURNS: C
  126.  */
  127. matrix4 *m4_add(matrix4 *A, matrix4 *B, matrix4 *C)
  128. {
  129.   word i;
  130.   register word j;
  131.  
  132.   for (i = 0; i < 4; i++) {
  133.     for (j = 0; j < 4; j++)
  134.       C->element[i][j] = A->element[i][j] + B->element[i][j];
  135.   }
  136.   return( C );
  137. }
  138.  
  139. /*----------------------------------------------------------------------------
  140.  * RETURNS: copy of A
  141.  */
  142. matrix4 *m4_copy(matrix4 *A)
  143. {
  144.   matrix4 *B = NEWTYPE(matrix4);
  145.   register word i, j;
  146.  
  147.   for (i = 0; i < 4; i++) {
  148.     for (j = 0; j < 4; j++)
  149.       B->element[i][j] = A->element[i][j];
  150.   }
  151.   return( B );
  152. }
  153.  
  154. /*----------------------------------------------------------------------------
  155.  * C = A * B
  156.  * RETURNS: C
  157.  */
  158. matrix4 *m4_mul(matrix4 *A, matrix4 *B, matrix4 *C)
  159. {
  160.   word i;
  161.   register word j, k;
  162.  
  163.   for (i = 0; i < 4; i++) {
  164.     for (j = 0; j < 4; j++) {
  165.       c->element[i][j] = 0;
  166.       for (k = 0; k < 4; k++)
  167.         C->element[i][j] += A->element[i][k] * B->element[k][j];
  168.     }
  169.   }
  170.   return( C );
  171. }
  172.  
  173. /*----------------------------------------------------------------------------
  174.  * RETURNS: new 4x4 matrix initialized to a[4][4]
  175.  */
  176. matrix4 *m4_new(double *a[4][4])
  177. {
  178.   matrix4 *A = NEWTYPE(matrix4);
  179.   register word i, j;
  180.  
  181.   for (i = 0; i < 4; i++) {
  182.     for (j = 0; j < 4; j++)
  183.       A->element[i][j] = *a[i][j];
  184.   }
  185.   return( A );
  186. }
  187.  
  188. /*----------------------------------------------------------------------------
  189.  * C = A - B
  190.  * RETURNS: C
  191.  */
  192. matrix4 *m4_sub(matrix4 *A, matrix4 *B, matrix4 *C)
  193. {
  194.   word i;
  195.   register word j;
  196.  
  197.   for (i = 0; i < 4; i++) {
  198.     for (j = 0; j < 4; j++)
  199.       C->element[i][j] = A->element[i][j] - B->element[i][j];
  200.   }
  201.   return( C );
  202. }
  203.  
  204. /*----------------------------------------------------------------------------
  205.  * B = transpose(A)
  206.  * RETURNS: B
  207.  */
  208. matrix4 *m4_transpose(matrix4 *A, matrix4 *B)
  209. {
  210.   register word i, j;
  211.  
  212.   for (i = 0; i < 4; i++) {
  213.     for (j = 0; j < 4; j++)
  214.       B->element[i][j] = A->element[j][i];
  215.   }
  216.   return( B );
  217. }
  218.  
  219. /*==========================================================================
  220.  *                                2D VECTORS
  221.  *--------------------------------------------------------------------------*/
  222.  
  223. /*----------------------------------------------------------------------------
  224.  * w = u + v
  225.  * RETURNS: w
  226.  */
  227. vector2 *v2_add(vector2 *u, vector2 *v, vector2 *w)
  228. {
  229.   w->x = u->x + v->x, w->y = u->y + v->y;
  230.   return( w );
  231. }
  232.         
  233. /*----------------------------------------------------------------------------
  234.  * w = c1(u) + c2(v)
  235.  * RETURNS: w
  236.  */
  237. vector2 *v2_combine(vector2 *u, double c1,
  238.                     vector2 *v, double c2, vector2 *w)
  239. {
  240.   w->x = (c1 * u->x) + (c2 * v->x);
  241.   w->y = (c1 * u->y) + (c2 * v->y);
  242.   return( w );
  243. }
  244.  
  245. /*----------------------------------------------------------------------------
  246.  * RETURNS: copy of u
  247.  */
  248. vector2 *v2_copy(vector2 *v)
  249. {
  250.   vector2 *w = NEWTYPE(vector2);
  251.  
  252.   w->x = v->x, w->y = v->y;
  253.   return( w );
  254. }
  255.         
  256. /*----------------------------------------------------------------------------
  257.  * w = u dot v
  258.  * RETURNS: w
  259.  */
  260. double v2_dot(vector2 *u, vector2 *v)
  261. {
  262.   return( (u->x * v->x) + (u->y * v->y) );
  263. }
  264.  
  265. /*----------------------------------------------------------------------------
  266.  * RETURNS: magnitude of v
  267.  */
  268. double v2_len(vector2 *v)
  269. {
  270.   return( sqrt(V2SquaredLength(v)) );
  271. }
  272.         
  273. /*----------------------------------------------------------------------------
  274.  * RETURNS: magnitude of v squared
  275.  */
  276. double v2_len_sqr(vector2 *v)
  277. {
  278.   return( (v->x * v->x) + (v->y * v->y) );
  279. }
  280.  
  281. /*----------------------------------------------------------------------------
  282.  * w = linear interpolation between lo and hi by amount alpha
  283.  * RETURNS: w
  284.  */
  285. vector2 *v2_lerp(vector2 *lo, vector2 *hi, double alpha, vector2 *w)
  286. {
  287.   w->x = LERP(alpha, lo->x, hi->x);
  288.   w->y = LERP(alpha, lo->y, hi->y);
  289.   return( w );
  290. }
  291.  
  292. /*----------------------------------------------------------------------------
  293.  * w = u * v (multiplication component-wise)
  294.  * RETURNS: w
  295.  */
  296. vector2 *v2_mul(vector2 *u, vector2 *v, vector2 *w)
  297. {
  298.   w->x = u->x * v->x;
  299.   w->y = u->y * v->y;
  300.   return( w );
  301. }
  302.  
  303. /*----------------------------------------------------------------------------
  304.  * RETURNS: transformed point p * projection matrix m
  305.  */
  306. point2 *v2_mul_by_proj(point2 *p, matrix3 *A)
  307. {
  308.   double w;
  309.   point2 q;
  310.  
  311.   q.x = (p->x * A->element[0][0]) +
  312.         (p->y * A->element[1][0]) + A->element[2][0];
  313.   q.y = (p->x * A->element[0][1]) +
  314.         (p->y * A->element[1][1]) + A->element[2][1];
  315.   w = (p->x * A->element[0][2]) +
  316.       (p->y * A->element[1][2]) + A->element[2][2];
  317.   if (w != 0.0)
  318.     q.x /= w, q.y /= w;
  319.   *p = q;
  320.   return( q );
  321. }
  322.  
  323. /*----------------------------------------------------------------------------
  324.  * v is negated
  325.  * RETURNS: v
  326.  */
  327. vector2 *v2_neg(vector2 *v)
  328. {
  329.   v->x = -v->x, v->y = -v->y;
  330.   return( v );
  331. }
  332.  
  333. /*----------------------------------------------------------------------------
  334.  * RETURNS: new 2D vector initialized to (x,y)
  335.  */
  336. vector2 *v2_new(double x, double y)
  337. {
  338.   vector2 *u = NEWTYPE(vector2);
  339.  
  340.   u->x = x, u->y = y;
  341.   return( u );
  342. }
  343.         
  344. /*----------------------------------------------------------------------------
  345.  * v is normalized (becomes unit length)
  346.  * RETURNS: v
  347.  */
  348. vector2 *v2_norm(vector2 *v)
  349. {
  350.   double len = V2Length(v);
  351.  
  352.   if (len != 0.0)
  353.     v->x /= len, v->y /= len;
  354.   return( v );
  355. }
  356.  
  357. /*----------------------------------------------------------------------------
  358.  * u = some orthogonal vector to v
  359.  * RETURNS: u
  360.  */
  361. vector2 *v2_ortho(vector2 *u, vector2 *v)
  362. {
  363.   v->x = -u->y;
  364.   v->y =  u->x;
  365.   return( v );
  366. }
  367.  
  368. /*----------------------------------------------------------------------------
  369.  * v is scaled to length newlen
  370.  * RETURNS: v
  371.  */
  372. vector2 *v2_scale(vector2 *v, double newlen)
  373. {
  374.   double len = V2Length(v);
  375.  
  376.   if (len != 0.0)
  377.     v->x *= newlen/len, v->y *= newlen/len;
  378.   return( v );
  379. }
  380.  
  381. /*----------------------------------------------------------------------------
  382.  * RETURNS: length of line segment pq
  383.  */
  384. double v2_seg_len(point2 *p, point2 *q)
  385. {
  386.   double dx = p->x - q->x;
  387.   double dy = p->y - q->y;
  388.  
  389.   return( sqrt(dx*dx + dy*dy) );
  390. }
  391.  
  392. /*----------------------------------------------------------------------------
  393.  * w = u - v
  394.  * RETURNS: w
  395.  */
  396. vector2 *v2_sub(vector2 *u, vector2 *v, vector2 *w)
  397. {
  398.   w->x = u->x - v->x, w->y = u->y - v->y;
  399.   return( w );
  400. }
  401.  
  402.  
  403. /*==========================================================================
  404.  *                                3D VECTORS
  405.  *--------------------------------------------------------------------------*/
  406.  
  407. /*----------------------------------------------------------------------------
  408.  * w = u + v
  409.  * RETURNS: w
  410.  */
  411. vector3 *v3_add(vector3 *u, vector3 *v, vector3 *w)
  412. {
  413.   w->x = u->x + v->x, w->y = u->y + v->y, w->z = u->z + v->z;
  414.   return( w );
  415. }
  416.  
  417. /*----------------------------------------------------------------------------
  418.  * w = c1(u) + c2(v)
  419.  * RETURNS: w
  420.  */
  421. vector3 *v3_combine(vector3 *u, double c1,
  422.                     vector3 *v, double c2, vector3 *w)
  423. {
  424.   w->x = (c1 * u->x) + (c2 * v->x);
  425.   w->y = (c1 * u->y) + (c2 * v->y);
  426.   w->z = (c1 * u->z) + (c2 * v->z);
  427.   return( w );
  428. }
  429.  
  430. /*----------------------------------------------------------------------------
  431.  * RETURNS: copy of u
  432.  */
  433. vector3 *v3_copy(vector3 *v)
  434. {
  435.   vector3 *w = NEWTYPE(vector3);
  436.  
  437.   w->x = v->x, w->y = v->y, w->z = v->z;
  438.   return( v );
  439. }
  440.  
  441. /*----------------------------------------------------------------------------
  442.  * w = u cross v
  443.  * RETURNS: w
  444.  */
  445. vector3 *v3_cross(vector3 *u, vector3 *v, vector3 *w)
  446. {
  447.   w->x = (u->y * v->z) - (u->z * v->y);
  448.   w->y = (u->z * v->x) - (u->x * v->z);
  449.   w->z = (u->x * v->y) - (u->y * v->x);
  450.   return( w );
  451. }
  452.  
  453. /*----------------------------------------------------------------------------
  454.  * w = u dot v
  455.  * RETURNS: w
  456.  */
  457. double v3_dot(vector3 *u, vector3 *v)
  458. {
  459.   return( (u->x * v->x) + (u->y * v->y) + (u->z * v->z) );
  460. }
  461.  
  462. /*----------------------------------------------------------------------------
  463.  * RETURNS: magnitude of v
  464.  */
  465. double v3_len(vector3 *v)
  466. {
  467.   return( sqrt(v3_len_sqr(v)) );
  468. }
  469.  
  470. /*----------------------------------------------------------------------------
  471.  * RETURNS: magnitude of v squared
  472.  */
  473. double v3_len_sqr(vector3 *v)
  474. {
  475.   return( (v->x * v->x) + (v->y * v->y) + (v->z * v->z) );
  476. }
  477.  
  478. /*----------------------------------------------------------------------------
  479.  * w = linear interpolation between lo and hi by amount alpha
  480.  * RETURNS: w
  481.  */
  482. vector3 *v3_lerp(vector3 *lo, vector3 *hi, double alpha, vector3 *w)
  483. {
  484.   w->x = LERP(alpha, lo->x, hi->x);
  485.   w->y = LERP(alpha, lo->y, hi->y);
  486.   w->z = LERP(alpha, lo->z, hi->z);
  487.   return( w );
  488. }
  489.  
  490. /*----------------------------------------------------------------------------
  491.  * w = u * v (multiplication component-wise)
  492.  * RETURNS: w
  493.  */
  494. vector3 *v3_mul(vector3 *u, vector3 *v, vector3 *w)
  495. {
  496.   w->x = u->x * v->x;
  497.   w->y = u->y * v->y;
  498.   w->z = u->z * v->z;
  499.   return( w );
  500. }
  501.  
  502. /*----------------------------------------------------------------------------
  503.  * q = p * A
  504.  * RETURNS: transformed point q
  505.  */
  506. point3 *v3_mul_by_matrix(point3 *p, matrix3 *A, point3 *q)
  507. {
  508.   q->x = (p->x * A->element[0][0])
  509.        + (p->y * A->element[1][0]) + (p->z * A->element[2][0]);
  510.   q->y = (p->x * A->element[0][1])
  511.        + (p->y * A->element[1][1]) + (p->z * A->element[2][1]);
  512.   q->z = (p->x * A->element[0][2])
  513.        + (p->y * A->element[1][2]) + (p->z * A->element[2][2]);
  514.   return( q );
  515. }
  516.  
  517. /*----------------------------------------------------------------------------
  518.  * RETURNS: transformed point p * projection matrix A
  519.  */
  520. point3 *v3_mul_by_proj(point3 *p, matrix4 *A)
  521. {
  522.   double w;
  523.   point3 q;
  524.  
  525.   q.x = (p->x * A->element[0][0]) + (p->y * A->element[1][0])
  526.       + (p->z * A->element[2][0]) + A->element[3][0];
  527.   q.y = (p->x * A->element[0][1]) + (p->y * A->element[1][1])
  528.       + (p->z * A->element[2][1]) + A->element[3][1];
  529.   q.z = (p->x * A->element[0][2]) + (p->y * A->element[1][2])
  530.       + (p->z * A->element[2][2]) + A->element[3][2];
  531.   w = (p->x * A->element[0][3]) + (p->y * A->element[1][3])
  532.     + (p->z * A->element[2][3]) + A->element[3][3];
  533.   if (w != 0.0)
  534.     q.x /= w, q.y /= w, q.z /= w;
  535.   *p = q;
  536.   return( q );
  537. }
  538.  
  539. /*----------------------------------------------------------------------------
  540.  * v is negated
  541.  * RETURNS: v
  542.  */
  543. vector3 *v3_neg(vector3 *v)
  544. {
  545.   v->x = -v->x, v->y = -v->y, v->z = -v->z;
  546.   return( v );
  547. }
  548.  
  549. /*----------------------------------------------------------------------------
  550.  * RETURNS: new 3D vector initialized to (x,y,z)
  551.  */
  552. vector3 *v3_new(double x, double y, double z)
  553. {
  554.   vector3 *v = NEWTYPE(vector3);
  555.  
  556.   v->x = x, v->y = y, v->z = z;
  557.   return( v );
  558. }
  559.  
  560. /*----------------------------------------------------------------------------
  561.  * v is normalized (becomes unit length)
  562.  * RETURNS: v
  563.  */
  564. vector3 *v3_norm(vector3 *v)
  565. {
  566.   double len = v3_len(v);
  567.  
  568.   if (len != 0.0)
  569.     v->x /= len, v->y /= len, v->z /= len;
  570.   return( v );
  571. }
  572.  
  573. /*----------------------------------------------------------------------------
  574.  * v is scaled to length newlen
  575.  * RETURNS: v
  576.  */
  577. vector3 *v3_scale(vector3 *v, double newlen)
  578. {
  579.   double len = v3Length(v);
  580.  
  581.   if (len != 0.0)
  582.     v->x *= newlen/len, v->y *= newlen/len, v->z *= newlen/len;
  583.   return( v );
  584. }
  585.         
  586. /*----------------------------------------------------------------------------
  587.  * RETURNS: length of line segment pq
  588.  */
  589. double v3_seg_len(point3 *p, point3 *q)
  590. {
  591.   double dx = p->x - q->x;
  592.   double dy = p->y - q->y;
  593.   double dz = p->z - q->z;
  594.  
  595.   return( sqrt(dx*dx + dy*dy + dz*dz) );
  596. }
  597.  
  598. /*----------------------------------------------------------------------------
  599.  * w = u - v
  600.  * RETURNS: w
  601.  */
  602. vector3 *v3_sub(vector3 *u, vector3 *v, vector3 *w)
  603. {
  604.   w->x = u->x - v->x;  w->y = u->y - v->y;  w->z = u->z - v->z;
  605.   return( w );
  606. }
  607.  
  608. /*==========================================================================
  609.  *                              USEFUL ROUTINES
  610.  *--------------------------------------------------------------------------*/
  611.  
  612. /*----------------------------------------------------------------------------
  613.  * binary greatest common divisor by Silver and Terzian.  See Knuth
  614.  * both inputs must be >= 0
  615.  */
  616. int gcd(int u, int v)
  617. {
  618.   int k, t, f;
  619.  
  620.   if (u < 0 || v < 0)
  621.     return( 1 );                       /* error if u < 0 or v < 0 */
  622.   k = 0, f = 1;
  623.   while ((u%2 == 0) && (v%2 == 0)) {
  624.     k++;
  625.     u >>= 1;
  626.     v >>= 1;
  627.     f *= 2;
  628.   }
  629.   if (u & 1) {
  630.     t = -v;
  631.     goto B4;
  632.   }
  633.   else
  634.     t = u;
  635. B3:
  636.   if (t > 0)
  637.     t >>= 1;
  638.   else
  639.     t = -((-t) >> 1);
  640. B4:
  641.   if (t%2 == 0)
  642.     goto B3;
  643.  
  644.   if (t > 0)
  645.     u = t;
  646.   else
  647.     v = -t;
  648.   if ((t = u-v) != 0)
  649.     goto B3;
  650.   return( u*f );
  651. }
  652.  
  653. /*----------------------------------------------------------------------------
  654.  * return roots of a*x^2 + b*x + c
  655.  * stable algebra derived from Numerical Recipes by Press et al
  656.  */
  657. int quadraticRoots(double a, double b, double c, double *roots)
  658. {
  659.   double d, q;
  660.   int count = 0;
  661.  
  662.   d = b*b - 4*a*c;
  663.   if (d < 0.0) {
  664.     *roots = *(roots+1) = 0.0;
  665.     return( 0 );
  666.   }
  667.   q = -0.5 * (b + SGN(b)*sqrt(d));
  668.   if (a != 0.0) {
  669.     *roots++ = q/a;
  670.     count++;
  671.   }
  672.   if (q != 0.0) {
  673.     *roots++ = c/q;
  674.     count++;
  675.   }
  676.   return( count );
  677. }
  678.  
  679.  
  680. /*----------------------------------------------------------------------------
  681.  * generic 1d regula-falsi step.  f is function to evaluate
  682.  * interval known to contain root is given in left, right
  683.  * returns new estimate
  684.  */
  685. double RegulaFalsi(double (*f)(), double left, double right)
  686. {
  687.   double d = (*f)(right) - (*f)(left);
  688.  
  689.   if (d != 0.0)
  690.     return( right - (*f)(right) * (right - left) / d);
  691.   return( (left+right) / 2.0 );
  692. }
  693.  
  694. /*----------------------------------------------------------------------------
  695.  * generic 1d Newton-Raphson step. f is function, df is derivative
  696.  * x is current best guess for root location. Returns new estimate
  697.  */
  698. double NewtonRaphson(double (*f)(), double (*df)(), double x)
  699. {
  700.   double d = (*df)(x);
  701.  
  702.   if (d != 0.0)
  703.     return( x - ((*f)(x) / d) );
  704.   return( x - 1.0 );
  705. }
  706.  
  707.  
  708. /*----------------------------------------------------------------------------
  709.  * hybrid 1d Newton-Raphson/Regula Falsi root finder
  710.  * input function f and its derivative df, an interval
  711.  * left, right known to contain the root, and an error tolerance
  712.  * Based on Blinn
  713.  */
  714. double findroot(double left, double right,
  715.          double tolerance, double (*f)(), double (*df)())
  716. {
  717.   double newx = left;
  718.  
  719.   while (ABS((*f)(newx)) > tolerance) {
  720.     newx = NewtonRaphson(f, df, newx);
  721.     if (newx < left || newx > right)
  722.       newx = RegulaFalsi(f, left, right);
  723.     if ((*f)(newx) * (*f)(left) <= 0.0)
  724.       right = newx;
  725.     else
  726.       left = newx;
  727.   }
  728.   return( newx );
  729. }
  730.  
  731. /*==============================  END-OF-FILE  =============================*/
  732.